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Abstract 

Neutron stars that are cold enough should have two or more superfluids/supercondutors in their 
inner crusts and cores. The implication of superfluidity/superconductivity for equilibrium and 
dynamical neutron star states is that each individual particle species that forms a condensate 
must have its own, independent number density current and equation of motion that determines 
that current. An important consequence of the quasiparticle nature of each condensate is the 
so-called entrainment effect, i.e. the momentum of a condensate is a linear combination of its own 
current and those of the other condensates. We present here the first fully relativistic modelling 
of slowly rotating superfluid neutron stars with entrainment that is accurate to the second-order 
in the rotation rates. The stars consist of superfluid neutrons, superconducting protons, and a 
highly degenerate, relativistic gas of electrons. We use a relativistic a - oj mean field model for the 
equation of state of the matter and the entrainment. We determine the effect of a relative rotation 
between the neutrons and protons on a star's total mass, shape, and Kepler, mass-shedding limit. 
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I. INTRODUCTION 



There is now an extended body of evidence in support of superfluidity in dense, nucleonic 
matter, such as that which is beheved to exist in neutron stars. Using conservative estimates 
one can argue for a Fermi temperature on the order of 10^^ K for neutrons in media with 
supranuclear densities, and that the transition temperature to a superfluid state is about 
10^ K. This is a signficiant fact since it is generally accepted that nascent neutron stars 
formed from supernovae should cool fairly quickly, and consequently their internal temper- 
atures should pass quickly through the transition value. Observational support for such 
tra.s,tions ,s suppl.ed by the well established ghtch phe„ome„o„ ■„ pulsars ^ These 
are rapid decreases in the rotational periods followed by a slow recovery |3], much too slow 
to be explained by ordinary fluid viscosity [J]. The best description is based on superfluid 
quantized vortices, and how they pin, unpin and then repin as the pulsar's rotation rate 
evolves We present here a fully relativistic formalism to model the rotational 

properties of superfluid neutron stars within a slow rotation approximation, which is an 
extension to second order in the rotation rates the previous work of Comer and Joynt ^. 

The bulk of the studies of neutron star superfluidity have been in the Newtonian regime, 
where the intent is not to be quantitatively descriptive but rather to gain qualitative under- 
standings. However, it is becoming increasingly apparent that general relativity is required 
to obtain even qualitative understandings. And even if the qualitative does not vary be- 
tween the two regimes, the number of examples is growing where general relativity can yield 
factors of two difference from Newtonian calculations (as opposed to "merely" 20% to 30% 
corrections). For example, recent modelling of supernovae has revealed that when gen- 
eral relativity is included the range of model parameters that exhibit multiple bounces is 
significantly smaller than the range found for the Newtonian case The use of general 
relativistic hydrodynamics has also led to predictions of the shock radius (during the shock 
reheating phase) being reduced by a factor of two and a corresponding increase by a factor 
of two for the inflow speed of the material behind the shock Certainly the need to use 

general relatrvity ,nust only be enhanced a= supernova renrnarrts beeonre nrore conrpact. 

Weber [11| provides an excellent overview of the many suggestions for the matter 
content of these remnants, the neutron stars. They range from the traditional neu- 
tron/proton/electron models to more exotic configurations with kaon or muon condensates 
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in the core, or hyperons, or even strange stars with absolutely stable u, d, s quark matter. In 
our study we will be somewhat conservative by considering typical traditional neutron star 
models, composed of superfluid neutrons, superconducting protons (with proton fractions 
on the order of 10%), and a highly degenerate relativistic gas of electrons. Each species 
extends from the center of the star all the way to the surface, although in principle one could 
consider a more realistic scenario wherein the protons and electrons extend out further than 



the neutrons 



13l ll^ and in that way mimic features of a crust. This same technique 
could be applied in neutron star cores, allowing for the possibility of alternating regions of 
ordinary fluid and superfluid. 

The local thermodynamic equilibrium of the matter will be modelled using a relativistic a 
- u mean field approach of the type that is attributed to Walecka and discussed in detail 
for neutron stars by Glendenning . We consider a relativistic approach to be important 
on two different levels. On the macroscopic level there is the need for general relativity 
that was discussed above. But on a microscopic level, recall that any fluid approximation 
for matter has built into it the notion of local, fluid elements. They are small enough that 
they can be considered to be points with respect to the rest of the star, and yet large enough 
to contain, say, an Avogadro's number of particles. At the densities expected for neutron 
stars, the local, Fermi levels for the nucleons can become high enough that the effective 
velocities of the nucleons with respect to their fluid elements are relativistic. As for the 

lieht, 



fluid elements themselves, they can also, in principle, approach speeds near that o: 
although in practice (e.g. for quasinormal mode [lJ,ll2| or slow rotation calculations 



they will typically have speeds a few percent of that of light. 

Because of the superfluidity of the neutrons, and superconductivity of the protons, the 
fluid formalism to be used differs fundamentally from the standard perfect fluid approach 
in that the neutrons can flow independently of the protons and electrons. There are thus 
two fluid degrees of freedom in the system, requiring two sets of fluid elements, one set for 
the neutrons and another for the charged constituents since the electromagnetic interaction 
very effectively "ties" the electrons to the protons. The matter description, therefore, must 
take into account two Fermi levels for the nucleons, and a displacement in momentum space 
between their respective Fermi spheres that will result when one fluid flows with respect 
to the other. The specification of a local thermodynamic equilibrium for the two fluids 
requires that the local neutron and "proton" (i.e. a conglomeration of the protons and 
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electrons) number densities be known as well as the local, relative velocity of the proton 
fluid elements, say, with respect to those of the neutrons. 

When the Fermi spheres for the nucleons are displaced with respect to each other there 
results an important effect for neutron star dynamics known as entrainment Sauls [2^ 

describes the entrainment effect using quasiparticle language, i.e. just because the neutrons 
are superfluid and protons superconducting does not mean they no longer feel the strong 
force. On the contrary, an individual neutron should be understood as being surrounded 
by a polarization cloud of other neutrons and protons. When this neutron moves, it will be 
accompanied by this cloud of nucleons. The net effect at the level of the fluid elements is 
that the momentum of a neutron fluid element is a linear combination of the neutron and 
proton number density currents. Parameters that are important for entrainment in neutron 



stars have been calculated in the Newtonian regime using a Fermi-liquid approach [21[ and 
in the relativistic regime via a a - u mean field model j^- Comer and Joynt ^ find that 
one of the key parameters that has been much used in superfluid neutron star modelling 
extends over a much larger range of values than what the Newtonian analysis of Borumand 
et al would imply 0|. 

Comer and Joynt also used their formalism to obtain first order rotational corrections to 
a superfluid neutron star's equilibrium state, which for general relativity means determining 
the frame-dragging, and the angular momentum. We will extend those calculations to 
the second order in the rotation rates, and thereby determine rotational corrections to the 
metric, distribution of particles, the total mass, shape, and the Kepler, mass-shedding limit. 
We construct sequences of equilibrium configurations where members of the sequence have 
the same relative rotation between the neutrons and protons but are distinguished by their 
central neutron number densities. We also construct sequences where the central neutron 
number density is fixed and the relative rotation is allowed to vary (cf . Q] for discussions 
on when to expect a relative rotation). One other straightforward, but important, extension 
here of the work of Comer and Joynt is a proof that the matter coefficients obtained by them 
are sufficient for our second order calculations. 

In order to have a reasonably self-contained document, and to define all the variables, we 
review in Sec. |n]the general relativistic superfluid formalism and its application to slowly 
rotating neutron stars. It is in this section that we prove the matter coefficients obtained 
by Comer and Joynt are all that is needed for the extension to second order. In Sec. IIIII 
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we discuss the highlights of the relativistic a - uj model and its mean field limit. We also 
determine the model's slow rotation limit. In Sec. IIVI we join the slow rotation formalism 
with the mean field model and produce numerical solutions. After reviewing the main 
results, the final, concluding section discusses applications beyond those considered here 
and points out where the formalism should be improved. For convenience we have restated 
in the appendix results of Comer and Joynt for the various matter coefficients that are 



required input for the field equations. We use "MTW' 
geometrical units. 



22| conventions throughout and 



II. GENERAL RELATIVISTIC SUPERFLUID FORMALISM AND SLOW ROTA- 
TION 

A. The full formalism 



We will use the 
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brmalism developed by Carter, Langlois, and their various collaborators 
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30|. The fundamental fiuid variables consist of the 



conserved neutron and proton number density currents, from which are formed the 
three scalars ra^ = —n^n'^, = —p^p^, and = —p^n'^. Given a master function 
— A(n^,p^,x^) (i.e. the superfiuid analog of the equation of state), then the stress-energy 
tensor is 



where 



is the generalized pressure and 



Ti: = m5>i + n^^li,+p^X. 



m = K- n^fip - pPxp 



/ij. = Briy + Apu , Xu = Aril, + Cp„ 



(1) 



(2) 



(3) 



are the chemical potential covectors which also function as the respective momenta for the 
fiuid elements. The A, B, and C coefficients are obtained from the master function via the 
partial derivatives 



A 



OA 



B 



, dA 



C 



^dA 

Qp2 



(4) 



5 



The fact that the neutron momentum /i^, say, is not simply proportional to its number 
density current n'^ is a result of entrainment between the neutrons and protons, which we 
see vanishes if the A coefficient is zero. 

Finally the equations for the neutrons and protons consist of the two conservation equa- 
tions 

vx = o , V^p^ = 0, (5) 

and the two Euler equations 

n^V[^fi,]=0 , p'^V[^XH = 0, (6) 



where the square braces mean antisymmetrization of the enclosed indices. Comer j30| and 



Prix et al 



121 discuss in some detail why the assumption of separate conservation for the 



two fluids should be reasonable for slow rotation and quasinormal mode calculations. 



B. The slow rotation expansions 



Andersson and Comer [18[ have adapted to the superfiuid case the ordinary fluid slow 
rotation scheme originally developed by Hartle j^^]- The configurations are assumed to be 
axisymmetric, asymptotically fiat, and stationary, with the metric taking the form 

g^Axf'dx'' = -(^N^ - sin^OK [N^] ^) dt^ + Vdr^ - 2sm^eKN'^dt d(f)+ 

K {de^ + sm^ed(j)^) . (7) 

The neutrons and protons are assumed to be rigidly rotating about the symmetry axis, with 
rates fin and flp, respectively, and unit four-velocities written as 

<= , + . <^ , + ^. (8) 



where is the Killing vector associated with the stationarity, and 0^ with the axisymmetry. 
The slow-rotation approximation assumes the rotation rates Q„ and Q„ are small in the 

nn 

sense that they should respect the inequalities (cf. |18L l31l|) 

^]>r^]Jorfi„^], « (|)'^, (9) 
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where the speed of hght c and Newton's constant G have been restored, and M and R are 
the mass and radius, respectively, of the non-rotating configuration. Because GM/c^R < 1, 
the inequalities naturally imply QnR and QpR << c. The slow- rotation scheme has 

been shown, for instance in Prix et al 12], to be a very good approximation for the fastest 
known pulsar, and starts to fail by about 15% to 20% for stars rotating at their Kepler limit. 
In the same manner as Hartle, the metric is expanded like 

N = e^(")/2(l + /i(r,^)) , 
V = e^^'^ {l + 2v{r,e)) , 
K = r\l + 2k{r,e)) , 

N''' = uj{r) , (10) 

where u is understood to be of 0{Qn,p) and h, v, and k of 0{fl'^p). For later convenience 
we will also introduce 

Ln = uj — and Lp = uj — Qp . (11) 

The expansion for the neutron and proton number densities n and p, respectively, are written 
as 

n = no{r){l + r]{r,e)) , p = Po(r) (1 + $(r, 0)) , (12) 

where the terms t] and $ are understood to be of C(f2^p) and we have introduced the 
convention that terms with an "o" subscript are either contributions from the non-rotating 
background or quantities that are evaluated on the non-rotating background (e.g. = rioPo 
etc.). One can show furthermore that the metric corrections h, v, and k can be decomposed 
into "/ = 0" and "Z = 2" terms using Legendre polynomials for the angular dependence, i.e. 

h = ho{r) + h2{r)P2{cose) , 
V = vo{r) + V2{r)P2{cos9) , 

k = k2{r)P2icos9) , (13) 
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where P2{cos6) = {Scos'^O — l)/2. As well we can write for the matter corrections 

V = r]o{r) + r]2{r)P2{cos9) , $ = <l>o(r) + <l>2(r)P2(cos6') . (14) 



Lastly, Andersson and Comer [Ig] have introduced a coordinate transformation r — >■ 
r + ^{r,6) that maps constant energy surfaces (i.e. the level surfaces of Aq) for the non- 
rotating background into the rotationally modified constant energy surfaces. The mapping 
^ is given as 

fJ'onoVO + XoPo^O + ^A^oPo i^n - ^pf = K^o (15) 

for I = and for / = 2 

fJ'onoV2 + XoPo*2 - IT^AoUoPo i^n " ^pf = Ao6 (16) 

where /Xo = BoUo + AoPo and Xo = C^Po + AoPo- 
C. The slow- rotation equations 

No different in number, but different in form from the ordinary fluid case considered 
by Hartle, the Einstein/superfluid field equations reduce to four sets: (i) the non-rotating 
background that determines z/. A, rio, andpo, (h) the linear order that calculates the "frame- 
dragging" u, (iii) the I = second order equations that yield ^o, ''70, '^'o, ^o, and Vq, and (iv) 
the Z = 2 set at second order that determines ^2, V2, ^2, h2, V2, and ^2. For convenience we 
list the equations here and recall where appropriate the free parameters of the system and 
how they are specified during a numerical integration. 



1. The 0(0^ „) or background equat 



ions 



The background equations are 

1 - 



A' = Svrre-^Ao , (17) 



= -- — — + Syrre^^o , (18) 



= ^SLp^+SoK + ^/^o'^', (19) 
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= C°|^p;+^°|^< + lxo^', (20) 

where the ^ol^, Bq\^, and Cq\^ coefficients are obtained from the master function and will be 
discussed in much greater detail below (cf . Sees. IllDl 11111 and the appendix) . Throughout 
a prime " ' " will denote differentiation with respect to r. 

Regularity at the origin implies that A(0), A'(0), i^'{0), n'^{0) and Po{0) all vanish. The 
radius R is obtained from the condition that the pressure vanish on the surface (i.e. \l/o(-R) = 
0) and the background mass is given by 

M = -An / Ao{rydr . (21) 
Jo 

The parameters no(0) and Po(0) for the matter on the background are not independent 
because chemical equilibrium is imposed between the neutrons, protons, and a highly de- 
generate gas of relativistic electrons (cf. the appendix). 



2. The 0{0,n,p) or frame- dragging equation 

For 0{Qn,p) there is only the frame- dragging uj{r), which is determined from 
i_ (^^4g-(A+.)/2^g' _ i6^e(A-.)/2 _ ^ Wne^'-'^y^XoPo {^^ - n,) . (22) 

The key difference with the frame-dragging equation for an ordinary fluid is the "source" 
term on the right-hand-side. Any solution must be such that it joins smoothly to the 
exterior vacuum solution. Continuity implies that 

2 7 

L„(i^) = -^]„ + — , (23) 

where J is the total angular momentum of the system. Using also that the derivative at 
the surface is to be continuous implies 

Z^i?) = - |LUi?) . (24) 

When the frame-dragging equation is solved numerically, one searches for a central value 
I/„(0) that allows Eq. to be satisfied. 

With a solution in hand, the neutron and proton angular momenta, J„ and Jp, respec- 
tively, can then be calculated using 

Jn = -— / drr^e(^-^)/2 + A^oPo iS^n - ^p) (25) 

•J Jo 



and 



57r 



3 ./o 



R 



drr^e(A-.)/2 



(26) 



The total angular momentum J likewise follows since J = Jn + Jp- 



3. The O(O^ p) equations 



One should first note that at this order there will exist both a / = and / = 2 set of 
equations. The 0{Q'^ p), Z = set consists of Eq. (fT^ as well as 

dA 



(27) 



Xo Xo ^e^Xo \ 



dA 

° dp 



+ rioPo 



OA 



(28) 



+ 



2 / r \ ' r'^ / ~ ^ ^ 
87rA;eo - ^ Hr^o + — ^ 1^1 



^2 



(29) 



2 2 / , 1 



\j^onoL^ + XoPoLp J 2^ 



(30) 



Here 7„ and 7p are integration constants, which are determined from data given at the center 
of the star: 



h (r\\ ^ I /^Q-^olo Xo^olo ^ ^ 



f^oXo 



r=0 



r=0 



(31) 
(32) 
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It is convenient to define the new variable mo = vvq/ exp(A). Because we must have 
(^o(O) = 0, then rjo{Q) and $o(0) are related via 



(/ionor/o)L=o + (XoPo$o)L=o = 



(33) 



Regularity at the center of the star also requires that mo(0) = 0, but /io(0) and 770 (0), say, 
are freely specified. 

As for the / = 2 equations, in addition to Eq. (fTBj) . we have 

OA 







Kino 



V2 H — '^'2 - 77-77— [Ao + rio — 



1^0 



dn 



3e' 



Xo 



V2 



dA 

dp 



(34) 



noPo 



dA 



(35) 



V2 + /?.2 



rrr 



XoPo (fin - fip) ( -^^n. + 



AqTIqPq (fiji fip) ] 



1 I/' 

= - {V2 + h2) - ik2 + h2)' --{h2- V2) 



(36) 
(37) 



2 6 2 
— r"'2 — ^/^a T 











- - ]V2^ 




--) 




r J 


r / 



4 r 
/co ^A;2 



"2 ^2' 



L'J + 



Stt 



'^'o - Ao) /i2 + — ( yUonoI^^ + XoPo^j 



7^ 2 



3e 



(38) 



In contrast to the / = case, the condition that ^2(0) = at the center of the star leads 
to ?72(0) = = $2(0). Also note that Eq. can be used to remove V2 in terms of the 
iunctions /i2 and Ln,p- It is furthermore convenient to work with the new variable k = h2 + k2 
. Regularity at the center of the star implies that h2{r) ~ cir^ and k{r) ~ C2r^ as r — 0, 



where the constants ci and C2 are related by 



ci + 27r ( ^o(O) - -Ao(0) ) C2 = 



(39) 
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We can take advantage of an overall scale invariance of the field equations to remove the 
explicit dependence on both fi„ and Qp by dividing through by Qp, the net result being 
that only the ratio Qn/^p needs to be specified in advance. As for the other parameters, 
they are chosen so that the solutions obtained join smoothly to a vacuum exterior. For 
/ = 0, after Qn/Qp and 770(0) are specified, this means a search over ho{0) is performed until 
a smooth solution is achieved. For / = 2, the situation is slightly more involved, requiring 
that homogeneous and then particular solutions be obtained for the set (/i2, k) (see [la] for 
more details). By searching over ci, say, and also trying different linear combinations of 
homogeneous and particular solutions eventually smoothness with the vacuum exterior can 
be achieved. 

Having obtained a complete solution, the rotationally induced change to the mass can be 
determined using 

5M={R- 2M) vo{R) + — . (40) 

Also the rotationally induced changes to the neutron and proton particle numbers are ob- 
tained, respectively, from 



and 



6N„ 



6N„ 





[A' 


2" 


2r2 




2^ 


h - 
r 


^0 + J 



^ drr^e^/^ + vq + 



2_ 



2r^~ 



(41) 



(42) 



Finally, the Kepler frequency is calculated (cf. j32] and [1^) using the slow-rotation 
form of 



Nv 



K 

where the metric variables are those of Eq. ((7j) and 

^3/2^;' 12KN' 



(43) 



NK' 



+ 



NK' 



+ 



NK' 



(44) 



is the orbital velocity according to a zero-angular momentum observer at the equator (where 
all quantities are to be evaluated). In the slow rotation approximation, the Kepler limit is 
the solution to 



JQp 



M ] SM ^ {R + 3M){3R-2M) ^3 




2M 



AR^M^ 



3 2^0 + 6 

4 i?2 



+ aA\nl (45) 
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where we have made the scahng with Qp exphcit by introducing J = JQp, etc., and have 
also defined 

3(i?3 - 2M3) , / 2M\ 3i?^ - 3R^M - 2R'^M'^ - 8RM^ + QM^ , , 

4M3 ^"H'"^J"' mM\R-2M) ■ 

In solving Eq. ^Jwe fix first the ratio fi„/fip. If |fi„/fip| < 1, then we insert = and 

solve for Qp, otherwise we set = and then again solve for flp. 

D. Analytic expansion for the local matter content 

Before proceeding to the next section that describes the mean field approach, one final 
observation needs to be made. We see in the equations above the appearance of not only 
the entrainment, via the A coefficient, but also the relative velocity-dependent part of the 
entrainment, through the term dA/dx^ . A priori one needs a formalism for the matter that 
will determine this term, if the effects of a relative rotation between the neutrons and protons 
are to be consistently incorporated. The mean field formalism presented in Sec. Illll can do 
just that. However, before embarking on a full-fiedged calculation, it is very worthwhile to 
consider next the implications of the slow rotation expansion for the master function. 

The term dA/dx^ would seem to imply, at least in principle, that we need to know the 
master function to 0{x^). Consider, then, regions of the matter (i.e. fiuid elements) that 
are local enough that their spacetime can be taken to be that of Minkowski, so that we can 
write out the x"^ term explicitly as 

where Vn and Vp are the neutron and proton, respectively, three-velocities in the local 
Minkowski frame. We see from this that when the individual three-velocities Vn,p satisfy 
Vn,p/c « 1, then x"^ ~ np to leading order in the ratios f„/c and Vp/c. 

With this as motivation we write an analytic expansion for the master function of the 
form 0, Q 

oo 

A(r^^p^x2) = ^A,(n^p2) [x^ - npY . (48) 
The A, A% etc. coefficients that appear in the field equations can now be written as 

I Ai [X" — np' * 

i=2 



A = -Ai -^i \i [x'^ - npY ^ , (49) 
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B = A Z^^[x -^P) > (50 

n on n n ^-^ on ^ ' 

C = ^ - - > — - np) , 51 

p op p p ^—^ op 

< = -l^-El^(-'--^^r > (52) 

OpOn ^ opOn ^ ' 
1=1 

1=1 

^t--'^-t^i-'-M' . (54) 

i=l 

^ = - - H - [i - {x' - npy-' , (55) 

^ = - ^ i [x^ - np] - [^ - 1] nAi^ (a;^ - np)'"' , (56) 

^ = -2A, - (z - 1) A, {x' - npy-' . (57) 



We conclude from this expansion that the master function coefficients Aq, Ai, and A2 
uniquely determine the background values for all the matter coefficients — i.e. Ao, B^, etc — 
that enter the field equations. In particular, if we need to know A2 then we need to know 
the master function to at least O(x^). However, the particular combinations that contain 
dA/dx"^ that enter the field equations are seen to reduce to 

. OA dA - dXi ^ /, dXA , o \i-i 

. dA dA , dXx ^ A (9AA , n 

That is, these combinations involve only the Ai coefficient when evaluated on the background. 
This, in fact, makes perfect sense since the slow rotation approximation is to 0(0^ p), which 
naturally corresponds to the term in the master function proportional to — np. 
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III. THE a - Lo RELATIVISTIC MEAN FIELD MODEL 



The master function to be used has been obtained using a relativistic a — uj mean field 
model of the type described by Glendenning The Lagrangian for this system is given 

by 

L = Li, + L, + L^ + Lint , (60) 

where 

U = ^{tl^^d^ - m)il} , (61) 
= -\d,ad^a-\mla\ (62) 
= -^ujf.^uj'"' - ^mluj^u^ , (63) 



Lint = gafTipip - g^u^il)'^^i) , (64) 

Here m is the baryon mass, is an 8-component spinor with the proton components as 
the top 4 and the neutron components as the bottom 4, the 7^ are the corresponding 8x8 
block diagonal Dirac matrices, and cj^jy = d^Ui, — d^uj^. The coupled set of field equations 
obtained from this Lagrangian are to be solved in each fiuid element of the neutron star. 
The main approximations of the mean field approach are to assume that the nucleons can 
be represented as plane-wave states and that all gradients of the a and fields can be 
ignored. The coupling constants and g^j and field masses m„ and are determined, for 
instance, from properties of nuclear matter at the nuclear saturation density. Fortunately, 
in what follows, we only need to give the ratios = (g^/m^)^ and = {gu/^LoY- 

Of course, the main consideration is to produce a master function that incorporates the 
entrainment effect. Consider again fiuid elements somewhere in the neutron star. The 
fermionic nature of the nucleons means that they are to be placed into the various energy 
levels (obtained from the mean field calculation) until their respective (local) Fermi spheres 
are filled. The Fermi spheres are best visualized in momentum space, as in Fig. ^ The 
entrainment is incorporated by displacing the center of the proton Fermi sphere from that 
of the neutron Fermi sphere, also illustrated in Fig. The neutron sphere is centered on 
the origin, and has a radius Displaced an amount K from the origin is the center of the 
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FIG. 1: The neutron and proton Fermi spheres drawn in momentum space. The neutron sphere 
is centered at the origin and has radius kn, whereas the proton sphere has radius kp and its center 
has been displaced an amount K from the origin. 

proton sphere, which has a radius kp. The Fermi sphere radii and displacement {kn,kp, K) 
are functions of the local neutron and proton number densities and the local, relative velocity 
of the protons, say, with respect to the neutrons. 
Introducing the definitions 



0n = fi'c^o; 



■p 



= (pn + K , 



(65) 



then the master function takes the form 



2cl 



1 



{ni^-ml)-?>{^>^^k,m) , 



(66) 



where 



127r2 



1 




dk, [kl - 2ml - 20^ - 3A;2 - 40„A;,) [kl + <pl + ml + 20„A;,) 



1/2 



+2 {% + (t)nf + ml 



2)3/2j + r dk, [(A;2 _ 2ml - 20j - ^kl - A<Ppk,) 



{kl + (Pl + ml + 2<Ppk,) + 2{% + (l)pf + m- 




(67) 
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= m — -—^m^ 



/Hp prCn j 

■/or) V 



hp 

1 



dK[{K + <Ppf + m'X'^\ , (68) 



= , (69) 

/ = ' (70) 



1 /'fen / 

= ^ / c^A;, (fc, + 0„) ( [kl + ml + <Pl + 20„A;,] - 

^ — hn 

[{k. + <Pn? + mlY^') , (71) 

= ^ f[ dkz {K + 0p) ( [fcj + + 0j + 20pA;,] - 

[(fc. + ^)Vm^]'/') , (72) 



0„ = -ci(n^+p^) . (73) 

As applied to this system, the slow-rotation approximation means that K is small with 
respect to kn,p- In particular, we should keep terms up to and including Oi^K"^). 

Implementation of the mean-field approach is hampered by the fact that the most nat- 
ural variables for it are the momenta (A;„, fcp, K), whereas the general relativistic superfiuid 
formalism relies on (n^,p^,x^). The two sets are related to each other via the kinematic 
relations 



n 



[n^'Y-in^f , = (p°)'- (p^)' , a;2 = ny-nV, (74) 



where n°, n^, and are given above in Eqs. (|69p. (j7(jp . (|7ip. and (|72p . In principle 
one would specify values for the set (n^,p^,x^), and then use Eqs. (fTOI). (f7T|) . and (f72|) 
in Eq. (f7^ to determine the set (fc„, fcp, i^). One would then solve Eq. (jUHj) for the Dirac 
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effective mass m*. In practice we rewrite the field equations so as to solve directly for 
kp, K) . 

Earlier we discussed the slow-rotation approximation as it applies to the master function, 
and determined that we need to know the Aq and Ai coefficients. Comer and Joynt |8| have 
shown that 



Ao 



^ (2m - m^\J (m - - ^ (m^kp [2kp + nie] ^Jkf+m^- 



m„ln 



kp + A^^p + 



(75) 



Ai 



2 V + 



fc« + -o^* I 



37r2 



A;2 + m,|„' 



+ 



A;2 + m,|^ 



(76) 



where m*!^ is the solution to the transcendental equation 



m — m* 



1 I 21 

- m 

2 '° 



^ I 21 

2 '° 



-fc„ + a/ fc„02 + mA 



(77) 



and we have added to Aq the contributions due to the electrons (me = m/1836). This is 
necessary to obtain a central proton fraction Po(0)/(po(0) + 7^o(0)) of about 0.1, which is 
considered typical for neutron stars. 

In order to solve for the background, we note that 

kl 



37r2 



Po 



kp 
3^ 



(78) 



That is, we replace everywhere the background neutron and proton number densities with 
their respective Fermi momenta. Once the Fermi momenta are specified, we can then deter- 
mine the background value for the Dirac effective mass, i.e. m*!^, from Eq. (f77j) . However, 
Comer and Joynt found that numerical solutions were easier to obtain by turning Eq. (f77|) 
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Model I 


Model II 


[gjm^f = 12.684 


{g^/m^f = 8.403 


{gu^/m^f = 7.148 


{gu^/m^f = 4.233 


1/(0) = -2.316408 


zv(0) = -2.288385 


A;n(0) = 2.8 fm"^ 


A;„(0) = 3.25 fm"^ 


x„(0) = 0.101 


xJQ) = 0.102 


M = 2.509 Mq 


M = 1.996 M0 


R = 11.696 km 


R = 9.432 km 


r?o(0) = 


??o(0) = 



TABLE I: The canonical background models used here and by Comer and Joynt 



into a differential equation using the identity 



dm J 



ftp, 



(79) 



and where k'^ and k'^ are to be obtained simultaneously using Eq. ()2()j) . The partial deriva- 
tives of the Dirac effective mass can be found in the appendix, along with the other matter 
coefficients. 



IV. NUMERICAL RESULTS 

In this section we present numerical solutions to the slow rotation equations. The two 
canonical background configurations_^(see Table HJ) that will be used repeatedly are the same 
as those used by Comer and Joynt 8] . They are such that the neutron and proton number 
densities vanish on the same surface. They are also such that chemical equilibrium has been 
imposed, which implies that the proton number density at the center of the star is no longer 
a free parameter. We do not apply chemical equilibrium for the rotating configurations 
since Andersson and Comer have shown it is not consistent with rigid rotation. The 
two sets of parameter values for the mean field model given in Table H] represent the two 
extremes discussed by Glendenning Q]. For all of the solutions we will set ?7o(0) = which 
implies also that $o(0) = 0. 

To begin we will give plots of the various field radial profiles for varying relative rotation 
rates of VLn/^p = 0.7, 1.0, 1.3. These fall within the range of rates that were considered ear- 
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FIG. 2: The metric functions ho{r) — the 3 lower curves — and /i2(r) — the 3 upper curves — vs r 
for Models I (left panel) and II (right panel) of Table ^ The relative rotation rates are varied 
(17„/0p = 0.7,1.0,1.3). 



lier by Comer and Joynt [8], who determined the frame-dragging and the angular momenta 
(cf. their Figs. 4 and 5). Picking up at the second order, we present in Fig. |21a plot of ho{r) 
vs r for several relative rotation rates and for Models I (left panel) and II (right panel) of 
TablelH Similarly in Fig. Elwe show a plot of mo{r) = rt'o(r)/ exp(A(r)) vs r and in Fig. |3] 
we have k{r) = h2{r) + k2{r) vs r. Fig.^lshows the radial profiles of ^o{r) and ^2{r)- Figs.lHl 
and [7| show the radial profiles of rio{r) and 772 (r) and ^o{r) and ^2{r), respectively. (We do 
not show plots of V2{r) since it can be removed in terms of h2{r) and Ln^p{r).) 

Andersson and Comer j3l used a simplified equation of state (i.e. a sum of polytropes) 
and did not include entrainment and thus a direct comparison between our solutions and 
theirs would not be expected to yield the same details. That being said, there are some 
qualitative similarities. For instance, the signs of the / = fields versus the / = 2 fields are 
in agreement, as are the number of zeroes that occur in the fields. The most pronounced 
qualitative differences are found in the functions ^o{f) and ^2{r), especially near the surface 
of the star. This could perhaps be related to our use of an explicit term for the electrons 
in the equation of state. 

In Fig. ISlwe have graphed the total mass for proton rotation rates that equal the fastest 
known pulsar (i.e. Qp = 3900 rad/s) versus the background central neutron number density. 
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FIG. 3: The metric function mo(r') = rfo(r)/exp(A(r)) vs r for Models I (left panel) and II (right 
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FIG. 4: The combination k{r) = /i2(r) + A;2(r) vs r for Models I (left panel) and II (right panel) of 
TableQ] The relative rotation rates are varied (r2„/ilp = 0.7, 1.0, 1.3). 



The left panel is for Model I of Table U and the right panel is for Model II. In both plots 
the relative rotation rate Qn/^p of the neutrons with respect to the protons is also varied. 
Generally speaking, stars that are spun up without changing the total baryon number move 
upward and to the left in the figure. We find that we recover the so-called "supramassive" 
configurations first discussed by Cook et al 3^. These are stars whose central densities 
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FIG. 5: The radial displacements £,o{r) — the 3 upper curves — and ^2(^) — the 3 lower curves — vs 
r for Models I (left panel) and II (right panel) of Tabled The relative rotation rates are varied 
(nn/Qp = 0.7,1.0,1.3). 

go beyond the maximum allowed for stable, non-rotating configurations and yet still remain 
on the stable (upward) branches of their particular mass vs central density curves. These 
configurations are stablized by their rotation. 

Fig. ini gives plots of the ratio of the polar (i.e. Rp = R + C,{R, 0)) to equatorial (i.e. Re = 
R + ^{R, 7r/2)) radii, versus the relative rotation rate Qn/^p for the canonical backgrounds 
given m Models I and II of Table HI and assuming a proton rotation rate of the fastest known 



pulsar. Up to and including O^Ql^p) we have 



Rp 

Rp 



1 + 



2 R 



(80) 



We see a general overall quadratic behaviour in the relative rotation, which is expected for 
the slow rotation expansion. Also expected is the limiting value of Rp/ Re = 1 as Qn/^p 0. 
Since the neutrons and protons have been constrained to vanish on the same surface, we 

14l | whereby the neutrons, say, can extend out 



do not find configurations like those of jl^ . 
further than the protons. 

Fig. ^1 is a graph of the Kepler (or mass-shedding) limit Qk versus the relative rotation 
rate for the canonical background configurations of Table IH Two competing factors are 
the amount of mass of the star that is in neutrons versus the relative rotation rate. For 
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FIG. 6: The neutron density corrections no{r)r]Q{r) — the 3 upper curves — and no{r)r]2{r) — the 3 
lower curves — vs r for Models I (left panel) and II (right panel) of Table The relative rotation 
rates are varied (0„/J7p = 0.7, 1.0, 1.3). 

instance, for f2„/f2p > 1 the Kepler limit is seen to be essentially flat. This can be easily 
understood as a result of the fact that the neutrons represent nearly 90% of the mass of the 
system, and so the Kepler limit changes very little as the relative rotation is increased. Such 
behavior has also been seen in the analytical solution of Prix et al , and the numerical 
results of Andersson and Comer flS^. The clear difference with the previous work is the 
general decrease in the Kepler limit as the relative rotation is decreased below one. Prix 
et al and Andersson and Comer found that the Kepler limit increased monotonically as the 
relative rotation was decreased. They explained this as a result of the fact that most of 
the mass is in the neutrons, and as the relative rotation is decreased, the Kepler limit is 
approaching the non-rotating value. As one can see in Fig. eventually a minimum is 
reached, and beyond that the Kepler limit then starts to increase as the relative rotation 
is decreased. It just does not exhibit the monotonic increase of the earlier studies. The 
main differences with the two earlier studies are that here we use the mean field formalism 
for the equation of state — the previous two studies used simplified polytropes — and the 
entrainment is both relativistic and also has a non-trivial radial profile (cf. Figs. 2 and 3 of 
0]), whereas the relevant entrainment parameter is kept constant in Prix et al and taken to 
be zero in Andersson and Comer. 
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FIG. 7: The proton density corrections Po('')^o('^) — the 3 uppermost curves at r/R = 1 — and 
Po(?')^2(?') — the 3 lowermost curves at r/R = 1 — vs r for Models I (left panel) and II (right panel) 
of Tabled The relative rotation rates are varied (J7„/ilp = 0.7, 1.0, 1.3). 

V. CONCLUSIONS 

The earlier work of Comer and Joynt laid the foundations for a fully relativistic 
approach to entrainment in superfluid neutron stars. They applied their formalism to first 
order in the rotational rates, and determined the impact of a relative rotation on the frame- 
dragging and the angular momenta of the two fluids. We have expanded on this work by 
extending the calculations to second order in the rotation rates, and have determined for the 
first time the maximum mass, shape, and Kepler limit using a fully relativistic formulation 
for entrainment and allowing a relative rotation between the two fluids. 

An important extension of our work will be to include isotopic spin terms in the master 
function ^l_6j. These are important for incorporating symmetry energy effects which tend 
to force baryonic systems to have as many protons as neutrons. Obviously nuclei tend to 
have equal numbers of neutrons and protons, and so this is why a symmetry energy term is 



12 



an important addition to the equation of state. In terms of neutron stars, Prix et al 
have shown clearly that symmetry energy impacts rotational equlibria (e.g. ellipticity, the 
Kepler Limit, and moment of inertia) as the relative rotation is varied. 

The next application we have in mind is to study quasinormal modes on slowly rotating 
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FIG. 8: The total mass as a function of central neutron number density for Models I (left panel) 
and 11 (right panel) of Tabled (using i7p = 3900 rad/s). The relative rotation rates are also varied, 
i.e. ^n/Vlp = 0,0.7,1.0,1.3. 



backgrounds. Andersson and Comer [3^ (see also [13]) have shown that a suitably advanced 
gravitational wave detector should be able to see gravitational waves emitted during Vela 
or Crab-like glitches. Such observations could provide potentially unique information on 
the supranuclear equation of state and the parameters that govern entrainment. Another 
important application will be to study further the recently discovered two-stream instability 



36|. It is the direct analog for superfluids of the instability of the same name known to 



exist in plasmas 
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38[ . Finally, our results should also find some application in studies of 
vortex structure in neutron stars. Link 3^ has suggested that free precession in neutron 
stars is not compatible with the protons behaving as a Type II superconductor. This should 
be explored in more detail by allowing for the effect of entrainment. 
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Appendix 

For convenience, we list here the various matter coefficients that are required as input in 
the field equations: 



B\ 



+ 



37r2 



i^n'^p 



kl + 



kl 



1,3 
2 l^p 



fO. 1.3 



kp + m,l 



37r2 



1.2 L.3 



kl + m^, 



+ 



lipl^n 



k^ + 



Sn^k^p kl + m^: 



+ 



(81) 



(82) 



26 



22000 



20000 - 



18000 



16000 



14000 



12000 



Model I 
Model 11 







0.5 



1 

Q. I Q. 



1.5 



FIG. 10: The Kepler mass-shedding Hmit as a function of the relative rotation rate Jln/fip for the 
canonical background configurations of Models I and II of Table ^ 



•^0 



^0 




^3 ^kl + ml 











1.5 










^ dkp 



kp + m,^, 1^ 
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^2 t J- m \ 



= c,2, + 



7]- A.p -r ^j^. 



kp ^kl + mi 



(83) 
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where 



dkn 



dkr. 



4 "^*lo^n 



^ 'kl + m* 



3m — 2 m* I c' 



m*L 77 



/c^ + m* 



+ 



A;^ + m* I 



3m-2m*|^ c; 



m* 



TT 



1,3 



kl + m* 



+ 



l3 



+ m*|„' 



^3 
ftp 



kl + mA 



(87) 



(88) 



The two functions /Iq and Xo are the background values for the two chemical potentials: 



m 



Xo 



{kl + kl) + Jkl 



m,* 



(89) 



Chemical equilibrium for the background means /Uq = Xo + \/kp + m^ must be imposed. 
Given a value for ^,^(0) then fcp(O) can be determined. Using these then the initial value for 
the Dirac effective mass can be obtained from Eq. (f77j) so that Eq. (f7^ can be integrated. 
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